function [W]= bciCSPEstimate (dataTS,label,method)

% W = megCSPEstimate(dataTS,label) transforms two dataset into a common component space
%
% Purpose: for a better classification result the covariance matrices of
%          two experimental conditions have to be passed as the input to
%          the function. The dataset has to have a
%          channels-by-time-points-by-trial format that CSP can cope with.
%
% INPUT:   timeseries datasets of the form mentioned above.
%
% OUTPUT:  a new matrix containing the common spatial components of which
%          some has to be selected and projected onto the trials
%
% REQUIREMENTS: The statistics toolbox

% 100215 SD wrote it
% 100302 SD Revision: according to Blankertz 2008 the simultanous
%                     diagonalization was changed (line no 37). Afterwards
%                     the new matrix is passed into a new order with
%                     respect to the variances of columns in EigVecs. As a
%                     result EigVecs is the output

if nargin<3,
    method=1; 
end
L=unique(label);
if length(L)~=2,
    error('Number of conditions doesn''t match two.');
end
c1 = find(label==L(1)); 
c2 = find(label==L(2));

% f�r jeden Trial wird die covarianzmatrix gebildet
Rx=zeros(size(dataTS,1),size(dataTS,1),length(c1));
for k = 1:length(c1),
%    dataTS(:,:,c1(k))=dataTS(:,:,c1(k))-repmat(mean(dataTS(:,:,c1(k)),2),[1,size(dataTS,2)]);
    Rx(:,:,k)  = (dataTS(:,:,c1(k) )*dataTS( :,:,c1(k) )')/trace( dataTS( :,:,c1(k) )* dataTS( :,:,c1(k) )' );
end;
Rx = mean(Rx, 3);
Ry=zeros(size(dataTS,1),size(dataTS,1),length(c2));
for k = 1:length(c2),
%    dataTS(:,:,c2(k))=dataTS(:,:,c2(k))-repmat(mean(dataTS(:,:,c2(k)),2),[1,size(dataTS,2)]);
    Ry(:,:,k)  = (dataTS(:,:,c2(k) )*dataTS( :,:,c2(k) )')/trace( dataTS( :,:,c2(k) )* dataTS( :,:,c2(k) )' );
end;
Ry = mean(Ry, 3);

if method == 1,  % Blankertz
    % aus der Summe beider Matrizen werden Eigenvektoren und Eigenwerte
    % bestimmt und nach absteigenden Eigenwerten geordnet
    [EigVecs,EigVals] = eig(Rx,Rx+Ry);
    % [EigVals,ind]     = sort(diag(EigVals),'descend');
    % EigVecs           = EigVecs(:,ind);
    % variances         = []; 
    % for filt = 1:length(EigVecs);
    %    variances(filt) = var(EigVecs(:,filt));
    % end
    % [variances, ind] = sort(variances);
    % W          = EigVecs(:,ind);
    W          = EigVecs;
else % Ramoser
    % Ramoser equation (2)
    Ccomposite=Rx+Ry;

    % Sort eigenvalues in descending order
    [Ucomposite,Lambdacomposite] = eig(Ccomposite);
    [Lambdacomposite,ind] = sort(diag(Lambdacomposite),'descend');
    Ucomposite = Ucomposite(:,ind);

    % Ramoser equation (3) - Whitening transform
    P=sqrt(inv(diag(Lambdacomposite)))*Ucomposite';

    % Ramoser equation (4)
    S{1}=P*Rx*P';
    S{2}=P*Ry*P';

    % Ramoser equation (5)
    [B,D] = eig(S{1},S{2}); % Simultanous diagonalization
			% Should be equivalent to [B,D]=eig(S{1});

    [D,ind] = sort(diag(D)); B = B(:,ind);

    W=(B'*P); % Projection matrix

    for i=1:length(ind), W(i,:)=W(i,:)./norm(W(i,:)); end

    A=pinv(W); % Common spatial patterns
    W=W';
end

if any(imag(W)~=0),
    warning('CSP filter includes complex numbers.');
end
 
 
 
 
 
    